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Abstract 

We  present  a  conceptual  discussion  of  a 
numerical  method  to  simulate  the  penetra¬ 
tion  of  a  ballistic  gel  by  a  projectile.  The 
method  draws  on  the  computational  advan¬ 
tages  presented  by  asynchronous  variational 
integrators  (AVIs)  and  an  immersed  bound¬ 
ary  method  to  result  in  an  easily  paral- 
lelizable,  efficient  and  adaptive  algorithm. 
The  paper  focuses  on  the  algorithmic  ideas 
involved-  at  the  expense  of  details  and  a  sim¬ 
plified  model. 

1  Introduction 

In  this  paper,  we  discuss  some  ideas  we  have 
for  simulating  the  mechanics  of  penetration  in 
a  soft-tissue  simulant.  Such  simulations  are 
of  immediate  interest  in  evaluating  firearms 
and  bullets  (see  fig.  1)  and  in  modelling  the 
complex  wounding  process.  Parameters  such 
as  the  depth  of  penetration  of  a  projectile, 
its  trajectory  or  the  extent  of  the  resulting 
wound  in  the  ballistic  gel(tissue  simulant)  can 
be  estimated  to  aid  in  optimizing  armor  ma¬ 
terials  and  designs. 

The  numerical  algorithm  that  we  propose 
consists  in  adopting  a  finite  element  dis¬ 
cretization  of  the  gel,  with  an  explicit  asyn¬ 
chronous  time  stepping  for  the  elements,  a 
contact  algorithm  to  handle  the  impact  of  the 
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Figure  1:  Snapshot  of  a  cavity 

left  in  a  ballistic  gel  by  a  bullet. 
Source:http://www.naaminis.com/pix/25gel02.jpg 

projectile  and  an  adaptive  remeshing  tech¬ 
nique.  AVIs  enable  each  element  in  the  finite 
element  mesh  to  be  advanced  independently 
(in  time)  with  a  time  step  befitting  the  local 
dynamics.  As  a  result,  small  time  steps  can 
be  adopted  for  elements  close  to  the  zone  of 
impact  or  in  regions  of  high  stresses  to  resolve 
events  occurring  at  small  time  scales. 

The  quality  of  the  mesh  used  for  the  gel  de¬ 
teriorates  as  the  projectile  penetrates-  a  re¬ 
sult  of  elements  with  large  strains  and  poor 
aspect  ratios.  It  is  therefore  necessary  to 
adopt  a  fresh  discretization  for  the  gel.  This 
is  achieved  with  an  adaptive  remeshing  al¬ 
gorithm  based  on  an  immersed  boundary 
method-  the  boundary  of  the  deformed  gel 
is  immersed  in  a  simple  discretization  from 
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which  a  new  mesh  for  the  gel,  adapted  to 
suit  the  local  dynamics,  is  extracted.  The 
displacement  and  velocity  fields  are  used  to 
dictate  the  element  sizes  in  the  new  mesh. 
These  fields  are  transferred  to  the  new  mesh 
and  the  simulation  is  set  up  to  continue. 

While  discussing  these  new  ideas  and  eval¬ 
uating  them  using  simple  examples,  we  make 
a  number  of  simplifying  assumptions.  While 
we  will  highlight  them  in  the  sections  that 
follow,  the  most  prominent  among  them  is 
that  we  model  the  phenomenon  as  an  im¬ 
pact  of  a  rigid  projectile  on  an  elastic  gel. 
Penetration  is  an  extremely  complex  phys¬ 
ical  phenomenon  heavily  dependent  on  the 
materials  involved  in  testing.  Experiments, 
for  instance  firing  high  speed  bullets  at  steel 
plates  of  different  thicknesses  (see  [2]),  reveal 
large  deformations,  high  strain  rates,  plastic 
flow,  fracture,  temperature  changes  as  well 
as  micro-structural  alterations  such  as  forma¬ 
tion  of  shear  bands.  It  is  not  hard  to  imagine 
a  transfer  of  mass  (besides  momentum)  from 
the  projectile  to  the  target  or  even  adhesion. 

Utility  of  numerical  simulations  such  as 
ours  depends  on  how  much  of  such  physics 
we  can  model.  But  there  are  significant  nu¬ 
merical  challenges  to  be  overcome.  Identi¬ 
fying  surfaces  of  contact,  avoiding  interpen¬ 
etration  of  colliding  bodies  or  handling  ge¬ 
ometries  with  sharp  corners  is  still  a  cumber¬ 
some  task  in  many  state-of-the-art  contact  al¬ 
gorithms.  There  is  still  some  way  to  go  be¬ 
fore  numerics  start  yielding  insight  about  the 
phenomenon  or  answering  specific  questions- 
such  as  the  influence  of  surface  roughness  in 
such  problems.  In  this  respect,  this  article 
should  be  viewed  as  a  progress  report  based 
on  exciting  ideas  we  have  to  address  some 
of  the  computational  issues  arising  in  simu¬ 
lations  of  penetration. 

In  the  following  sections,  we  discuss  the 
numerical  algorithms  we  adopt  to  simulate 
a  simplified  penetration  problem-  an  asyn¬ 
chronous  variational  time  integrator,  a  con- 


fa)  Snapshot  of  a  simulation  of  the  impact 
of  an  ‘L’  shaped  beam  against  a  rigid  wall 
(not  shown)  using  AVI  and  the  contact  algo¬ 
rithm  described  next.  Elements  are  colored 
according  to  the  time  step  adopted,  with  blue 
representing  an  order  of  magnitude  smaller 
than  red.  The  computational  efficiency  of 
AVIs  enables  the  resolution  of  high  frequency 
modes  seen  in  such  problems 


Time 

(b)  Energy  conservation  characteristics  of  AVIs 
demonstrated  in  the  impact  problem.  The  energy 
is  conserved  nearly  exactly,  as  the  inset  shows, 
with  only  a  minor  drop  in  value  through  impact, 
which  can  be  regulated  through  the  selection  of 
the  time  step  for  the  elements  in  contact. 

Figure  2:  AVI  with  contact  for  dynamics  simu¬ 
lations. 
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(a)  An  idealization  of  the  proposed  contact  algo¬ 
rithm.  The  contact  constraint  is  given  by  the  con¬ 
dition  g  <  0.  When  contact  occurs  (g  =  0)  with 
the  rigid  projectile,  the  normal  component  of  the 
momentum  of  the  gel  is  reversed  while  keeping  the 
tangential  component  unaltered. 


(b)  Interpenetration  is  permitted  in  this  algorithm. 
Even  though  contact  may  be  detected  after  the  gel 
has  penetrated  the  projectile,  this  can  be  made  neg¬ 
ligible  by  adopting  small  time  steps  for  elements 
close  to  the  contact  region. 


Figure  3:  Illustration  of  the  contact  algorithm. 


tact  algorithm  and  an  adaptive  remeshing 
strategy.  We  will  focus  on  the  main  ideas 
rather  delving  into  specifics,  which  can  be 
found  in  the  references  provided. 

2  Asynchronous  Varia¬ 
tional  Time-Integrators 

AVI  possesses  two  properties  most  desirable 
in  dynamic  simulations, 

•  being  variational,  it  has  outstanding  mo¬ 
mentum  and  energy  conservation  prop¬ 
erties,  and 

•  being  asynchronous,  it  has  the  distin¬ 
guishing  feature  of  permitting  an  inde¬ 
pendent  selection  of  time  steps  for  each 
element  in  a  mesh. 

Fig. 2  shows  a  snapshot  of  the  simulation  of 
the  impact  of  an  ‘L’  shaped  beam  against  a 
rigid  wall  using  AVI  and  the  almost  exact  en¬ 
ergy  conservation  of  the  system. 

With  traditional  time  integrators,  the  time 
step  for  the  entire  mesh  is  dictated  by  the 


minimum  over  the  range  of  stable  time  steps 
for  elements  in  the  mesh.  This  can  be  partic¬ 
ularly  hindering  in  high  velocity  contact  sim¬ 
ulations.  Local  stiffening  of  the  material  and 
adaptively  refined  meshes  can  cripple  long¬ 
time  simulations.  With  AVI,  each  element’s 
time  step  need  only  satisfy  its  own  stabil¬ 
ity  criterion,  dependent  on  the  local  sound 
speed,  material  velocity  and  element  size, 
rather  than  the  most  stringent  one  for  the 
mesh.  See  [6]  for  a  detailed  derivation  and 
analysis  of  the  algorithm  and  [3]  for  imple¬ 
mentation  details. 

3  Contact  Algorithm 

Next,  we  discuss  a  contact  algorithm  formu¬ 
lated  to  handle  the  impact  of  the  projectile 
on  the  gel.  The  projectile  is  approximated  as 
a  rigid  body  and  the  gel  as  a  soft  hyperelas¬ 
tic  material.  The  rigid  body  approximation 
is  adopted  mainly  for  the  significant  algorith¬ 
mic  simplifications  it  yields,  though  this  can 
be  loosely  justified  by  the  contrast  in  proper¬ 
ties  of  the  projectile  and  the  gel. 
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(a)  Partition  of  the  domain  for  parallel  pro¬ 
cessing. 


(b)  An  elastic  (Neohookean)  block  impacted  by 
a  rigid  sphere. 


(c)  A  slice  through  the  center  of  the  block  show¬ 
ing  the  large  deformation.  The  color  contours 
correspond  to  the  magnitude  of  velocity,  de¬ 
creasing  in  magnitude  from  blue  to  red. 

Figure  4:  A  snapshot  of  the  impact  of  a  rigid 
sphere  on  an  elastic  block  modeled  using  the  con¬ 
tact  algorithm  described. 


Traditionally,  contact  has  been  modeled  as 
a  constrained  optimization  problem  in  which 
an  energy  functional  is  rendered  stationary 
subject  to  the  constraint  of  no  interpene¬ 
tration.  In  the  context  of  the  current  (dis¬ 
cretized)  problem,  such  a  constraint  can  be 
expressed  in  the  form  g(x)  <  0,  where  x  rep¬ 
resents  the  set  of  nodal  positions  of  the  gel 
and  g  is  a  function  that  is  strictly  positive 
whenever  a  node  lies  in  the  interior  of  the  pro¬ 
jectile,  and  strictly  negative  when  all  nodes 
lie  in  the  exterior.  Note  that  g  is  implicitly 
a  function  of  time  as  well  via  the  positions  of 
the  projectile  and  the  gel. 

With  this  description,  collision  occurs 
when  g(x)  =  C  >  0.  Upon  detecting  a  col¬ 
lision  (at  the  end  of  a  time  step  of  some  el¬ 
ement),  the  component  of  the  momentum  of 
the  gel  normal  to  the  surface  g(x)  —  C  is 
reversed  while  leaving  the  tangential  compo¬ 
nent  unchanged,  see  fig.  3.  Note  that  no  pro¬ 
jection  onto  the  boundary  g(x)  =  0  is  per¬ 
formed  to  remove  any  penetration  that  has 
occurred,  although  this  is  one  possibility.  In 
this  way,  the  gel  effectively  sees  the  contact 
“wall”  wherever  it  falls  at  the  end  of  a  time 
step,  instead  of  at  a  fixed  location.  An  unde¬ 
sirable  consequence  of  permitting  such  inter¬ 
penetration  is  the  dissipation  of  energy.  How¬ 
ever,  the  time  steps  for  elements  approach¬ 
ing  contact  with  the  projectile  can  be  taken 
to  be  very  small  (leaving  the  time  steps  else¬ 
where  unaltered)  so  that  interpenetration  and 
consequently  the  energy  loss  is  small.  The 
almost  exact  energy  conservation  in  the  ex¬ 
ample  of  the  L-shaped  beam  shown  above 
demonstrates  precisely  this.  Another  exam¬ 
ple  of  a  rigid  sphere  impacting  a  block  is 
shown  in  fig.  4. 

In  comparison  with  the  above  algorithm, 
an  approach  of  computing  the  exact  time  of 
contact  for  each  node  close  to  the  projec¬ 
tile  would  prove  prohibitively  expensive.  A 
penalty  approach  is  commonly  used  to  min¬ 
imize  interpenetration,  whereby  the  poten- 
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tial  in  the  (interior  of  the)  projectile  is  made 
large,  approaching  that  of  a  rigid  body,  so 
that  it  is  energetically  unfavorable  for  the 
gel  to  cross  the  boundary  of  the  projectile. 
While  this  has  the  advantage  that  an  inequal¬ 
ity  constraint  need  not  be  explicitly  handled, 
determining  how  stiff  the  penalty  potential 
should  be  while  maintaining  accuracy  and 
avoiding  ill-conditioning  issues  is  often  non 
trivial.  Likewise,  Lagrange  multipliers  have 
also  been  used  to  impose  the  contact  con¬ 
straint  at  the  cost  of  adding  to  the  number 
of  unknowns  in  the  system,  and  making  it 
implicit.  A  comprehensive  review  of  contact 
algorithms  can  be  found  in  [9,  10]. 

4  A  Preliminary  Penetra¬ 
tion  Algorithm  With 
Adaptive  Remeshing 

In  principle,  the  explicit  time  integration  and 
contact  algorithms  are  sufficient  to  simulate 
the  long-term  behavior  of  the  dynamic  sys¬ 
tem  under  consideration.  However,  as  the 
simulation  progresses,  elements  in  the  gel  get 
significantly  deformed  resulting  in  poor  as¬ 
pect  ratios,  as  shown  in  fig.  5.  This  neces¬ 
sitates  a  smaller  time  step  to  ensure  the  sta¬ 
bility  of  the  integrator  and  finally  results  in 
an  unfeasible  time  step.  The  remedy  lies  in 
adopting  a  better  discretization  for  the  cur¬ 
rent  configuration.  This  is  precisely  what  we 
discuss  next  with  the  aid  of  a  2 D  example. 
Most  of  the  ideas  can  be  extended  to  the  3D 
case  in  a  straightforward  manner. 

The  adaptive  remeshing  algorithm  broadly 
consists  of  the  following  steps: 

(a)  determining  when  the  quality  of  the 
mesh  of  the  gel  has  deteriorated  enough 
to  warrant  a  new  discretization, 

(b)  computing  an  indicator  of  how  fine  the 
mesh  need  locally  be  to  represent  all  cur¬ 


ia)  Mesh  undergoing  large  deformation. 


(b)  Close-up  view  of  the  elements  in  the  region 
of  contact. 

Figure  5:  A  2D  example  of  a  rigid  ball  impacting 
a  square  in  which  poorly  shaped  elements  result 
from  large  deformations  of  the  gel  impacted  by  a 
rigid  ball.  This  renders  explicit  time  integration 
algorithms  unstable  unless  extremely  small  time 
steps  are  used. 

rent  velocity,  strain  stress  and  displace¬ 
ment  fields, 

(c)  creating  a  new  mesh  for  the  current  con¬ 
figuration  of  the  mesh  with  element  sizes 
roughly  as  specified  by  the  indicator  and 

(d)  transferring  fields  (such  as  displacements 
and  velocities)  via  interpolation  to  ones 
on  the  new  mesh. 

No  doubt,  each  of  the  above  steps  is  a  dif¬ 
ficult  question  and  there  have  been  numer¬ 
ous  solutions  proposed,  each  tailored  to  suit 
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(a)  We  adopt  an  ad-hoc  local 
refinement  indicator  based  on 
identifying  regions  of  high  cur¬ 
vature  in  the  displacement  or 
velocity  fields.  Shown  here  is 
the  curvature  (in  degrees)  for 
the  vertical  component  of  the 
piecewise  linear  velocity  field. 
Note  that  change  in  angles  can 
be  as  large  as  90°. 


(d)  The  gel  is  considered  to  be 
immersed  in  the  triangulation 
of  the  square.  To  extract  a 
mesh  for  the  gel,  the  linear  in- 
terpolant  of  the  signed  distance 
function  to  the  boundary  of  the 
gel  is  computed  over  all  nodes 
of  the  triangulation. 


(b)  The  remeshing  strategy 
involves  the  construction  of 
a  quadtree  over  a  square  do¬ 
main  encompassing  the  gel. 


(e)  A  discretization  for 
the  gel  is  determined  by 
trimming  (and  subdividing 
quadrilateral  cuts)  elements 
with  the  zero  level-set  of  the 
linearized  signed  distance 
function. 


(c)  A  Delaunay  triangulation  is 
built  over  the  convex  hull  of  the 
centers  of  the  quadtree  nodes. 
Note  that  this  triangulation  is 
sufficiently  refined  in  regions 
where  the  indicator  is  large  in 
the  domain  of  the  gel. 


(f)  Displacement  and  velocity 
fields  are  transferred  to  the  new 
mesh  via  interpolation.  Shown 
here  is  the  vertical  component 
of  the  velocity  field 


Figure  6:  Pictorial  description  of  the  adaptive  remeshing  algorithm.  The  distinguishing  feature  of 
this  algorithm  is  its  simplicity-  it  avoids  many  of  the  cumbersome  operations  (like  edge  flipping, 
mesh  smoothing)  commonly  associated  with  meshing.  It  extends  to  3 D  in  a  rather  straightforward 
manner  and  can  be  automated  quite  easily. 
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Figure  7:  Deformation  of  a  square  block,  impacted  by  a  circular  ball,  computed  using  the  remeshing 
strategy  described.  The  block  was  remeshed  four  times  (after  fixed  intervals)  without  considering 
adaptivity.  Shown  also  are  the  contours  of  the  vertical  component  of  the  velocity. 


a  particular  application.  Shown  in  fig.  6  is  a 
pictorial  description  of  the  adaptive  remesh¬ 
ing  algorithm  that  we  propose.  The  highlight 
is  the  remeshing  strategy  that  draws  from 
a  discontinuous  Galerkin  based  immersed 
boundary  method  [5,  8]  we  recently  formu¬ 
lated  for  solid  mechanics.  The  idea  is  to 
immerse  the  boundary  of  the  domain  to  be 
meshed  (the  gel  in  this  case)  in  a  simple  dis¬ 
cretization  and  extract  a  mesh  that  approxi¬ 
mately  fits  the  boundary.  In  this  way,  espe¬ 
cially  in  3D,  the  difficulties  associated  with 
automatic  mesh  generation  are  quite  easily 
overcome.  The  rationale  behind  permitting 
an  approximate  mesh  is  that  it  would  be  fu¬ 
tile  to  attempt  a  perfectly  conforming  dis¬ 
cretization  considering  the  numerical  approx¬ 
imations  involved  and  that  the  geometries 
are  already  only  approximately  represented. 
Moreover,  the  approximations  in  the  domain 
are  very  small-  decreasing  quadratically  with 
the  new  mesh  size  along  the  boundary  of  the 
gel.  It  is  important  to  note  that  even  though 
with  this  approach  elements  near  the  bound¬ 
ary  are  often  smaller  than  those  in  the  inte¬ 
rior,  thanks  to  AVI,  they  run  with  their  own 
time  step,  without  affecting  the  overall  time 
step  of  the  system. 


the  robustness  of  the  procedure,  the  Achilles 
heel  of  every  existing  such  approach. 

5  Summary 

Penetration  problems  are  of  great  practical  as 
well  as  academic  interest,  playing  host  to  a 
whole  range  of  physical  phenomena.  Numer¬ 
ical  methods  are  a  key  tool  to  understanding 
the  underlying  physics.  But  there  are  signif¬ 
icant  challenges  to  be  overcome  to  gain  real¬ 
istic  insights  (rather  than  being  misled  by  al¬ 
gorithmic  artifacts).  The  numerical  method 
presented  here  is  a  step  in  this  direction.  Be¬ 
sides  its  computational  efficiency  and  capabil¬ 
ity  to  simulate  long-term  dynamics,  we  also 
try  to  thoroughly  understand  the  approxima¬ 
tions  involved.  With  the  (afore  mentioned) 
numerical  tools,  modern  computational  re¬ 
sources  and  more  realistic  models  of  the  pene¬ 
tration  phenomenon,  simulations  such  as  ours 
can  be  expected  to  reveal  exciting  details. 


The  simulations  shown  above  and  in  fig.  7 
are  only  preliminary  ideas.  These  have  to  still 
be  refined  and  precisely  analyzed  to  enhance 
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